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Abstract: The statistics of recurrence times in broad areas have been 
reported to obey universal scaling laws, both for single homogeneous regions 
(Corral, 2003) and when averaged over multiple regions (Bak et al.,2002). 
These unified scaling laws are characterized by intermediate power law asymp- 
totics. On the other hand, Molchan (2005) has presented a mathematical 
proof that, if such a universal law exists, it is necessarily an exponential, in 
obvious contradiction with the data. First, we generalize Molchan's argument 
to show that an approximate unified law can be found which is compatible 
with the empirical observations when incorporating the impact of the Omori 
law of earthquake triggering. We then develop the full theory of the statis- 
tics of inter-event times in the framework of the ETAS model of triggered 
seismicity and show that the empirical observations can be fully explained. 
Our theoretical expression fits well the empirical statistics over the whole 
range of recurrence times, accounting for different regimes by using only the 
physics of triggering quantified by Omori's law. The description of the statis- 
tics of recurrence times over multiple regions requires an additional subtle 
statistical derivation that maps the fractal geometry of earthquake epicenters 
onto the distribution of the average seismic rates in multiple regions. This 
yields a prediction in excellent agreement with the empirical data for rea- 
sonable values of the fractal dimension d ^ 1.8, the average clustering ratio 
n ~ 0.9, and the productivity exponent a ~ 0.9 times the 6-value of the 
Gutenberg-Richter law. Our predictions are remarkably robust with respect 
to the magnitude threshold used to select observable events. These results 
extend previous works which have shown that much of the empirical phe- 
nomenology of seismicity can be explained by carefully taking into account 
the physics of triggering between earthquakes. 
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1 Introduction 



The concept of "recurrence time" (also called recurrence interval, inter- 
event time or return time) is widely used for hazard assessment in seismology. 
The average recurrence time of an earthquake is usually defined as the num- 
ber of years between occurrences of an earthquake of a given magnitude 
in a particular area. Once chosen a probabilistic model for the distribu- 
tion of recurrence times (often a Poisson distribution in earlier implementa- 
tions, evolving progressively to more realistic distributions), an estimation 
of the probability for damaging earthquakes in a given area over a specific 
time horizon can be derived. For instance, the Working Group on Califo- 
nia Earthquake Probabilities (2003) assessed a 0.62 probabihty of a major 
damaging earthquake striking the greater San Francisco Bay Region over the 
following 30 years (2002-2031). Such calculations are based on classifications 
of known active faults and on assumptions about the organization of seis- 
micity on these faults, guided by geological, paleoseismic (see for instance 
Sieh, 1981), geodetic evidence (Global Earthquake Satellite System, 2003) as 
well as earthquake patterns extracted from recent seismic catalogs. Starting 
with the simple assumption that faults are independent and carry their own 
characteristic earthquake (Schwartz and Coppersmith, 1984), seismologists 
and other scientists working on seismic hazards are realizing that more com- 
plex models are needed to take into account the interaction, coupling and 
competition between faults (Lee et al., 1999) which translates into a rich 
phenomenology in the space-time-magnitude organization of earthquakes. 

A new view is now emerging that recurrence times should be considered 
for broad areas, rather than for individual faults, and could provide impor- 
tant insights in the physical mechanisms of earthquakes. For instance, Wyss 
et al. (2000) have proposed to relate an indirect measurement of local re- 
currence times within faults to geometrical asperities and stress field, thus 
significantly broadening the standard notion of return time. Inspired by the 
scaling approach developed in the physics of critical phenomena (Sornette, 
2004), Kossobokov and Mazhkenov (1988), Bak et al. (2002) and Christensen 
et al. (2002) have proposed a unified scaling law combining the Gutenberg- 
Richter law, the Omori law and the fractal distribution of epicenters to de- 
scribe the distributions of inter-event times between successive earthquakes 
in a hierarchy of spatial domain sizes and magnitudes in Southern Califor- 
nia. Corral (2003, 2004a,b, 2005a) has refined and extended these analyses 
to many different regions of the world and has proposed the existence of a 
universal scaling law for the probability density function (PDF) H{t) of re- 
currence times (or inter-event times) r between earthquakes in a given region 
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S: 



H{t) ^ A/(Ar) 



(1) 



The remarkable proposition is that the function f{x), which exhibit differ- 
ent power law regimes with cross-overs, is found almost the same for many 
different seismic regions, suggesting universality. The specificity of a given 
region seems to be captured solely by the average rate A of observable events 
in that region, which fixes the only relevant characteristic time 1/A for the 
recurrence times. 

The interpretation proposed by Bak et al. (2002), Christensen et al. 
(2002) and Corral (2003, 2004a,b, 2005a) is that the scaling law (P) reveals 
a complex spatio-temporal organization of seismicity, which can be viewed 
as an intermittent flow of energy released within a self-organized (critical?) 
system, for which concepts and tools from the theory of critical phenomena 
can be applied (Corral, 2005b). This view has been challenged by Lindman 
et al. (2005) who stressed several methodological caveats (see Corral and 
Christensen (2006) for a reply). Livina et al. (2005) have additionally noticed 
that the marginal (or mono-variate) PDF of inter-event times gives only a 
partial description of the time sequence of earthquakes in a given region, as 
it is found to be a function of preceding inter-event times. There is thus a 
memory between successive earthquakes which influences the distribution of 
the inter-event times, a conclusion which is well-known to most seismologists. 

It is fair to state that there is at present no theoretical understanding 
of these empirical results and in particular of expression The situation 
becomes even more interesting with the recent mathematical demonstration 
by Molchan (2005) that, under very weak and general conditions, the only 
possible form for f{x), if universality holds, is the exponential function, in 
strong disagreement with the observations reported by Bak et al. (2002) 
and Corral (2003, 2004a,b, 2005a). In addition, from a re-analysis of the 
seismicity of Southern California, Molchan and Kronrod (2005) have shown 
that the unified scaling law (P) is incompatible with multifractality which 
seems to offer a better description of the data. 

Here, our purpose is to show how all the above can be simply understood 
and reconciled from the standard known statistical laws of seismicity: 

• (i) the Gutenberg- Richter distribution ~ l/E-^~^^ (with /? ^ (2/3)6 ^ 
2/3) of earthquake energies E (Knopoff et al., 1982); 

• (ii) the Omori law ~ l/t^ (with p ^ 1 for large earthquakes) of the rate 
of aftershocks as a function of time t since a mainshock (Utsu et al., 
1995); 
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• (iii) the productivity law ~ E"' (with a ~ 2/3) giving the number 
of earthquakes triggered by an event of energy E (Helmstetter et aL, 
2005); 

• (iv) the fractal (and even probably multifractal (Ouillon et al., 1996)) 
structure of fault networks (Davy et al, 1990) and of the set of earth- 
quake epicenters (Kagan and Knopoff, 1980). 

The question we address can be summarized as follows: are the statistics on 
inter-event times described in (Bak et al, 2002; Corral, 2003, 2004a,b, 2005a; 
Livina et al.. 2005) really new in the sense that they reveal some information 
which is not contained in the above laws (i-iv), as claimed by these authors? 
Or can they be derived from the known statistical properties of seismicity, 
so that they are only different ways of presenting the same information? 

In order to address this question, we use the simplest possible model which 
combines these four laws, framed as a stochastic point process of past earth- 
quakes triggering future earthquakes: the Omori law (ii) is taken to describe 
the conditional rate of activation of new earthquakes, given all past earth- 
quakes; the Gutenberg-Richter law (i) describes the distribution according 
to which the magnitude of a new earthquake is independently determined; 
the productivity law (iii) gives the weight of the contribution of a given past 
earthquake in the production of new earthquakes. In order to take into ac- 
count the fractal geometry of earthquake catalogs, the new earthquakes can 
be positioned on a fractal geometry. The elements (i-iii) actually constitute 
the basic building block of the benchmark model of seismicity, known as the 
Epidemic- Type Aftershock Sequence (ETAS) model of triggered seismicity 
(Ogata, 1988; Kagan and Knopoff, 1981), whose main statistical properties 
are reviewed in (Helmstetter and Sornette, 2002). Various versions have been 
or are currently applied by different groups to describe observed and forecast 
future seismicity (Console et al., 2002; 2003a,b; 2006; Console and Murru, 
2001; Gerstenberger et al, 2005; Reasenberg, R A. and Jones, 1989; 1994; 
Steacy et al., 2005). 

The common characteristics of this class of models is to treat all earth- 
quakes on the same footing such that one does not assume any distinction 
between foreshocks, mainshocks and aftershocks: each earthquake is con- 
sidered to be capable of triggering other earthquakes according to the three 
basic laws (i-iii) mentioned above. This hypothesis is based on the realization 
that there arc no observable differences between foreshocks, mainshocks and 
aftershocks (Jones et al., 1999; Helmstetter and Sornette, 2003a; Helmstetter 
et al., 2003), notwithstanding their usual classification based in retrospect 
analysis of realized seismic sequences: Earthquakes preceding the mainshock 
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are called foreshocks and earthquakes following the mainshock are called af- 
tershocks but no physical differences between these earthquakes are known. 
For instance, a mainshock (classified as the largest event in a given space- 
time window) may be re-classificd as a foreshock if a larger event later follows 
it; and what would have been an aftershock becomes a mainshock if larger 
than all its preceding earthquakes. 

In this paper, we use specifically the ETAS model in the version proposed 
by Ogata (1988). The ETAS model assumes that earthquakes magnitudes are 
mutually statistically independent and drawn from the Gutenbcrg-Richter 
(GR) probability Q{m). The GR law gives the probability that magnitudes 
of triggered events are larger than a given level m (the relationship between 
magnitude m and energy is m oc (2/3) luio E). We shall use the GR law in 
the form 

Q(m2) = g(mi)10-'('"^-"^i) , (2) 

which emphasizes its scale invariance property. Here 6 ~ 1 and mi, m2 
are arbitrary magnitudes. We parameterize the (bare) Omori law (Sornette 
and Sornette, 1999; Hclmstetter and Sornette, 2002) for the rate of triggered 
events of first-generation from a given earthquake as 

with 9 >0. One may interpret $(t) as the PDF of random times of indepen- 
dently occurring first-generation aftershocks, triggered by some mainshock 
which happened at the origin t = 0. The last ingredient is the productivity 
law, which we write in the following convenient form for further analysis, 

p{m) = K 10°("^-"*°) , (4) 

where the factor k, will be related below to physically observable quantities 
such as the average branching ratio n (which is the average number of trig- 
gered earthquakes of first generation per triggering event). The magnitude 
mo is a cut-off introduced to regularize the theory (Helmstetter and Sornette, 
2002). It can be interpreted as the smallest possible magnitude for earth- 
quakes to be able to trigger other earthquakes (Sornette and Werner, 2005a). 
Several authors have shown that the ETAS model provides a good description 
of many of the regularities of seismicity (Console et al., 2002; 2003a,b; 2006; 
Console and Murru, 2001; Helmstetter and Sornette, 2003a,b; Helmstetter et 
al., 2005; Gerstenberger et al., 2005; Ogata, 1988; 2005; Ogata and Zhuang, 
2006; Reasenberg, P. A. and Jones, 1989; 1994; Saichev and Sornette, 2005; 
2006a; Steacy et al., 2005; Zhuang et al., 2002; 2004; 2005). 



5 



Our main result is that, according to Occam's razor, the previously men- 
tioned results on universal scaling laws of inter-event times do not reveal 
more information than what is already captured by the well-known laws 
(i-iii) of seismicity (Gutenberg- Richter, Omori, essentially), together with 
the assumption that all earthquakes are similar (no distinction between fore- 
shocks, mainshocks and aftershocks), which is the key ingredient of the ETAS 
model. Our theory is able to account quantitatively for the empirical power 
laws found by Bak et al. and Corral, showing that they result from subtle 
cross-overs rather than being genuine asymptotic scaling laws. We also show 
that universality does not strictly hold. 

The organization of the paper is the following. In section 2, we discuss 
in detail Molchan's derivation (Molchan, 2005) that, if the scaling law (Q) 
holds true, then it should be exponential in the sense that f{x) = e~^. We 
extend Molchan's argument by developing a simple semi- quantitative theory, 
showing that the presence of the Omori law destroys the self-similarity of the 
PDF of recurrence times. Nevertheless, if the exponent of the Omori law Q 
is close to 1, i.e., ^ -C 1, then, the PDF of recurrence times approximately 
obeys a non-exponential scaling law which can fit well the empirical data. 
Section 3 extends further the discussion of section 2 by proposing simplified 
models of aftershock triggering process, in the goal of demonstrating the di- 
rect relation between the Omori law and the corresponding scaling law for 
the PDF of the recurrence times for an arbitrary region. These discussions 
allow us to stress the generality of the curve that we propose to the dis- 
agreement between Molchan's result and empirical data. They also prepare 
us to understand better the full derivation using the technology of generat- 
ing probability functions (GPF) applied to the ETAS model. In section 4, 
we analyze the ETAS model of triggered seismicity with the formalism of 
GPF and establish the main general results useful for the following. As a 
check for the formalism, we obtain the statistical description for the number 
of observable events which occur within a given region and during the time 
window [t,t + t]. Section 5 exploit the general results of section 4 to obtain 
predictions on the PDF of recurrence times between observable events in a 
single homogeneous region characterized by a well-defined average seismic 
rate. A short account of some of these results is presented in (Saichev and 
Sornette, 2006c). Section 6 combines the previous results on single regions 
to obtain predictions of the PDF of inter-event times between earthquakes 
averaged over many different regions. In a first part, we construct different 
ad hoc models of the distribution of the average seismic rates of different 
regions. In a second part, we propose a procedure converting the fractal ge- 
ometry of epicenters into a specific distribution of the average seismic rates 
of different regions. The knowledge of this distribution allows us to compute 
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the PDF of inter-event times averaged over many regions which is compared 
with Bak et al. (2002) 's and Corral (2004a) 's empirical analyses. We find 
a very good agreement between our prediction and the empirical PDF of 
inter-event times. 

2 Generalized Molchan's relation 

We start by addressing the puzzle raised by the demonstration by Molchan 
(2005) based on a probabilistic reasoning that, under very weak and general 
conditions, the only possible form for f{x) in (^, if universality holds, is 
the exponential function, in contradiction with the strongly non-exponential 
form of the reported seemingly universal unified PDF for inter-event 
times. Here, we briefiy reproduce Molchan's argumentation and then gen- 
eralize it to show that the presence of the Omori law, while destroying the 
exact unified law, gives nevertheless an approximate unified law fitting rather 
precisely the real data. 

Let -P(t) be the probability that there are no observable events during 
the time interval [t, t + t] within some region S. Let us call A the average 
rate of observable events within the region S. Let us assume the existence of 
a unified law 

P(r) = (/.(Ar) , (5) 

valid for any region. In expression (0), (p{x) is the universal scaling function, 
which is assumed to be the same for all regions. Following Molchan's argu- 
ment, if an unified law ^ holds true for any region, then it should be valid 
for the region made of the union of two disjoint regions characterized respec- 
tively by the rates Ai and A2. If the triggering processes of earthquake in 
those two regions are statistically independent, then the following functional 
equation should be true 

¥'((Ai + A2)r) =y.(Air)¥;(A2r) . (6) 

Posing -^(x) = \iinp{x), the equation (0) is equivalent to 

^^"^ = ^(Arib^")+^(A;i^") • 

It is well-known that the solution of equation (jZj) defines the class of linear 
funtions tplx) = —ax for arbitrary a's, leading to the exponential form 

^(x) = e-'^- . (8) 
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For the PDF (jSJ, there is a normahzation condition that allows us to fix 
the constant a. We need first to recall a few general results of the theory of 
point processes (Daley and Vere- Jones, 1988) Let the random times 

• • • t-i-,to-, ti, . . .tk . ■ ■ tk+i > tk > tk-i ... (9) 

form a stationary point process with an average rate equal to A. Then, 
the PDF H{t) of the inter-event times tk+i — tk can be obtained from the 
following relation 

where -P(t) is the already mentioned probability of absence of events within 
the interval + r]. The normalization condition for -ff(r) reads 



A dr 



T = 



=1. (11) 

x=0 ^ ' 



Substituting with (jH)) yields a = 1. Thus, following Molchan's reasoning, if 
a unified law (Q) for the PDF of recurrence times exists, then it should be 
exponential 

/(x) = ^'\x) = e-^ . (12) 

We have already mentioned that the law (fT^ derived by Molchan (2005) 
contradicts the observations (Bak et al., 2002; Corral, 2003, 2004a,b, 2005a). 
Let us thus generalize Molchan's reasoning, by assuming that the probability 
P(r) depends, not only on the dimensionless combination Ar, but addition- 
ally on the average rate A. This means that we assume that the probability 
P(r) can expressed as 

P(r) = v9(Ar,A) , (13) 
for some universal function ip. With this assumption, equation (jH)) is replaced 

by 

y?((Ai + A2)r, Al + A2) = ^{XiT, Ai)^(A2r, A2) , (14) 
and expression (jZj) is changed into 

ij{x, Al + A2) = (x;^^' ^1) + ^ ( aTTa;''' ' ^^^^ 

where again ip = Imp and x is a dimensionless variable. It is easy to check 
that the solution of (fT3j) has the form 

^(x. A) = -X g (^j^ , (16) 
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where g{y) is a function which is arbitrary except that it should be such 
that ip{x, A) is monotonically decreasing with respect to x. Expression (fTBj) 
translates into 

P(r) = e-^"3(") . (17) 

It is clear from expression ()17|1 that g{T) should be dimensionless, hence 
its argument must be dimensionless. Since g is independent of A, the only 
possibility to make g{T) dimensionless is that it depends on another time 
scale c such that g can be written. 

9{r) = go(^^) ' ^^^^ 

where go{y) is a function with dimensionless argument. This implies that 
expression (fTTjl has actually the form 

P(r)=exp(-Ar(7o(^)) , (19) 

which generalizes ((HI). Note that the later is recovered for the special case 
where the function go is a constant (which has then to be unity by the 
condition of normalization discussed above). 

The relevance of a time scale c in addition to the inverse rate 1/A is 
actually part of Omori's law. Consider the often used pure power form of 
the Omori law 

$(t) = k r^-^ . (20) 

Then, necessarily, for 6* 7^ 0, a time scale c is needed so that $ has the 
dimension of the inverse of time, i.e., 

k = koc^ (21) 

where fco is a dimensionless constant. 

It appears natural on physical grounds (and will be justified in our cal- 
culations below) that the Omori law has an influence on the form of the 
probability P{t). This suggests that the physical origin of the time scale 
c (and hence of the deviation (fT^ from Molchan's law (jH))) lies in Omori's 
law. In other words, this argues for the fact that the function g^ in (fT^ 
actually derives from Omori's law. We thus obtain the key insight that the 
existence of the Omori law implies the absence of a universal scaling law for 
the PDF of inter-event times! 

In our calculations presented below using the ETAS model, we will show 
that the Omori law gives the following specific structure for the function 

goiy) = a + hy-r (22) 
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Substituting it into (fT^ and using the dimensionless variable x = Ar, we 
obtain 

P(r) = ip{x, e) = exp (-ax - e^/ix^"^) , e = Ac . (23) 

A remarkable property of the distribution (f{x, e) is that it changes very 
slowly for 6 <^1, even if the average rate A changes by factors of thousands. 
This characteristic property of ip{x, e) is at the origin, as we show below, of 
the essentially non-exponential approximate unified law for the PDF (Q) of 
recurrence times, which provides excellent fit to the PDF's of real data. 

3 Simplified model of the impact of Omori's 
law on the PDF of inter-event times 

In the preceding section, we have suggested that the Omori law may be at 
the origin of both (i) a breaking of the unified PDF of the inter-event times 
from the exponential function derived by Molchan and (ii) an approximate 
universal law different from the exponential. In the following sections, we will 
present a rigorous analysis of the PDF of recurrence times in the framework 
of the ETAS model, which confirms this claim. In the present section, we 
provide what we believe is a simple intuitive understanding of the roots of the 
observed approximate unified law, by using a simplified model of recurrence 
time statistics which takes into account the impact of the Omori law. 

We consider a synthetic Earth in which two types of earthquakes oc- 
cur. Spontaneous events are triggered by the driven tectonic forces. In turn, 
these spontaneous events (also called "immigrants" in the jargon of branching 
processes) may trigger their "aftershocks" (called more generally "triggered 
events"). We put quotation marks around "aftershocks" to stress the fact 
that those "aftershocks" may be larger than their mother event. They are 
not necessarily the aftershocks of the nomenclature of standard seismology. 
These triggered events of first generation may themselves be the sources of 
triggered events of the second generation and so on. We denote u the Pois- 
son rate of the observable spontaneous events. The average number of first- 
generation "aftershocks," triggered by some spontaneous event, defines the 
key parameter n of the theory, often called the branching parameter. For the 
process to be stationary and not explode, we consider the sub-critical regime 
n < 1 (Helmstetter and Sornette, 2002). Note that n is by definition the aver- 
age number of first-generation events triggered by any arbitrary earthquake. 
It has also the physical meaning of being the fraction of triggered events in 
a catalogue including both spontaneous sources and triggered events over all 
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possible generations (Helmstetter and Sornette, 2003c). Given the average 
rate u) of the spontaneous sources, the average rate of all events, including 
spontaneous and all their children over all generations, is 07(1 + n + + ...), 
that is, 

A = -^. (24) 

\ — n 

In reality, one does not have the luxury of observing all earthquakes. Cat- 
alogues are complete only above a minimum magnitude which depends on 
the density of spatial coverage of the network of seismic stations and on 
their quality. One thus typically observes only a small subset of all occurring 
earthquakes, with the vast majority (of small earthquakes) being hidden from 
detection. Here, we formulate a simplifying hypothesis, justified later by a 
full rigorous calculation with the ETAS model, that, due to the statistical 
independence of the magnitudes of earthquakes, the ratio of the number of 
observable spontaneous events to the number of all observable events is ap- 
proximately independent of the magnitude threshold m. Assuming this to be 
true, then the relation ()24|) also holds between the rate ujim) of observable 
spontaneous events and the rate A(m) of the total set of observable events, 
whose magnitudes all exceed the threshold level m. 

We take into account the influence of both the spontaneous sources and of 
the triggered events on the probability -P(t) of abscence of events within time 
interval [t,t -l- r], by postulating the following generalized Poisson statistics 

P(r) ^exp(-cjr-cjA(r)) . (25) 

The first term lot in the exponential is nothing but the standard Poisson 
contribution that no spontaneous events fall inside the time interval [t, t-|-r]. 
The second term co'A(r) describes the contribution resulting from all the 
events which have been triggered before t. Neglecting for simplicity the 
difference of productivity of different events, we use the simplified form 

^A(r) = \{t')dt' , A(t') ^ *(^' - , (26) 

where ^{t) is the normalized Omori law described by relation (jH)). We assume 
that the Omori law is self-averaging, so that we can replace A(t') by its 
average value 

(A) (t' -t) = r $(t' - t")dt" = -^a(t' - t) . (27) 

1 — n J-00 1 — n 

Here, ujn/{l — n) is the average number per unit time of triggered events 
occurring before the time window [t, t + t] and 



a{t) = I mdt' = j--^ ■ (28) 
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Replacing in (j2SI) A(t') by (A)(f - t) yields 



A(r) 



n 



1 — n 



Air) 



where 



A{t) 



a{t)dt = 

1 — 6' 



T 

1 + - 

C 



In particular, if r 3> c, we have 



A( 



un c / T 
1-nl-e \c 



t:^ c 



(29) 



(30) 



(31) 



Thus, the simplified model ()25|) for the probability that no events occur in 
the time interval [t, t + t] takes the form 



P(r) = exp —LOT 



ujn 
1 — n 



A{r] 



which, for t ^ c, reads 

F(r) f« exp —ojT 



Lon c / T 



t:^ c . 



l-nl- 

Using relation ()24|1 and introducing the dimensionless parameters 

X = Ar , e = Ac , 

we obtain finally 



P(r) ^ (p{x, e) = exp —(1 — n)x 



ne 

1-e 



X 



1-6 



(32) 



(33) 



(34) 



(35) 



This expression fl35|) coincides with the previously proposed form ()23p with 
the correspondence a = 1 — n and h = n/(l — 9). 

By substituting into (fTn|l . we obtain that H{t) is described by the 
relation (^, where 



f{x) = ne'Ox-'-' + 



1 — n + ne x' 



Lp{x,i 



(36) 



The PDF /(x) (jHUj) depends on the scale L of the region under ob- 
servation, via the dependence of the parameter e on the rate A which itself 
depends on L. Similarly, via the same chain of dependence, f{x) is also a 
function of the threshold magnitude m of observed events. However, the 
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PDF f{x) given by ()36|) is very slowly varying with the average rate A for 
9 1, which is typical of aftershock sequences. When A varies by large 
factors up to thousands, changes insignificantly. Indeed, let us denote by 
Ai and A2 the average rates corresponding to two magnitude threshold levels 
TUi and 7712- Using the GR law, the parameter changes by the factor p 
times = pe^) where 



Consider for instance 777,2 — fTT-i = 4, for which the average rate changes by 
the huge factor 10~^ while the p ~ 0.76 still remains rather close to 1 for 
e = 0.03. Fig. □ plots the PDF f{x) given by §^ for n = 0.9, 6 = 0.03 
and for different magnitude thresholds m. This figure demonstrates the slow 
dependence of f{x) with respect to magnitude threshold level m. 

While the plots of the approximate scaling law ()3(ij) shown in Fig. ^ differ 
significantly from Molchan's exponential unified law, they are close to the 
parameterization of the empirical PDF proposed by Corral (2004a) 



with the parameters 

7 = 0.67 ±0.05, 5 =1.05 ±0.05, a = 1.64 ±0.15, (39) 

and with a normalization factor C close to 1. Fig. El shows a typical PDF 
f{x) obtained from our expression together with Corral's fitting curve 
()38|1 . illustrating their closeness for not too small x's. The two expressions 
depart for very small x < 0.01, for which our simplified model distribution 
f{x) given by (jHUIl exhibits a power asymptotic f{x) ~ x^^^^, which is a 
direct consequence of the Omori law. This asymptotic term results explicitely 
from the first term in the factor multiplying (p{x, e) in expression ()36|1 . This 
asymptotic is absent in Corral's fitting function fc{x) given by (jHEl), although 
it actually exists in real data plots (Corral, 2004a). We will emphasize this 
point later when we discuss the full theory and its comparison with the 
empirical data. 

4 Description of the statistics of observable 
events in the context of the ETAS model 

In the preceding sections, we have made plausible and intuitively appeal- 
ing the possibility that an approximate unified law indeed exists for the PDF 




(37) 




(38) 
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of the inter-earthquake times, by using a simphfied transparent model of the 
statistics of recurrence times. In this section, we develop an accurate math- 
ematical description of the statistical properties of observable events, in the 
framework of the ETAS model presented in the introduction. The technology 
developed here will then be used in the next section to quantify precisely the 
prediction of the ETAS model of the PDF of inter-event times and to derive 
different approximations. 

Our analysis is based on our previous calculations (Saichev and Sornette, 
2006a), which showed that, for large domain sizes (L ~ tens of kilometers or 
more), one may neglect the impact of aftershocks triggered by events that 
occurred outside the considered spatial domain S, while only considering the 
events within the space domain which are triggered by sources also within 
the same domain. In this approximation validated by precise calculations 
in (Saichev and Sornette, 2006a), the analysis of the statistical properties 
of time sequence of earthquakes can be reduced to the study of a time-only 
version of the ETAS model, in which the dependence on the spatial extension 
of a given area can be taken into account via the dependence on L of the 
rate uj of spontaneous observable events. 

As the ETAS model has an exact representation in terms of branching 
clusters (Hawkes and Oakes, 1974; Helmstetter and Sornette, 2002), the 
main mathematical tool to analyze the ETAS model in this context is the 
Generating Probability Function (GPF) of the random number R[t, r, m) of 
observable events within the time interval [t, t + r] 

(](^;r,m) = (z^^*'"'™)) . (40) 

Here and below, the angle brackets (...) represent the procedure of statistical 
averaging. We assume that the process of earthquake triggering is stationary, 
so that the GPF given by ()40|) does not depend on t but only on the duration 
r of the time interval [t,t + r]. Let P{r,T,m) be the probability that the 
number R{t, r, m) of observable events is equal to r. Then, one may rewrite 
the definition (pn|) in the equivalent form 

oo 

^l{z;T,m) = ^P{r;T,m)z'^ . (41) 

r=0 

It follows that the probability that no event occurs in a time interval of 
duration r is equal to 

P{t) = P{r = 0,T,m) = n{0;T,m) . (42) 

To determine the GPF Q{z; r, mo) associated with the total number of 
events inside the time window [t,t + r], we proceed as follows. We partition 
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the t-axis in small intervals Ik = [tk,tk + A] as shown in Fig. Eland then 
group them in three categories: the first one I^^ G {—oo,t) lies before the 
time window [t, t+r] of interest, the second one G [t, t+r] groups all small 
intervals inside the time window [t, t + r] and the third one J^g G (t + r, oo) 
lies after the time window [t, t + r]. Let us now determine in turn the impacts 
of the spontaneous events occurring within each of these three intervals on 
the GPF of the number of events inside the time window [t,t + r]. Let 
Qki{z; , mo) be the GPF of the number of spontaneous events occurring in 
the first class /^^ of small intervals. It is equal to 

oo 

ek,{z;mo) = Y.nr,^)z\ (43) 

r=0 

where V{r, A) is the Poissonian probability that r spontaneous events occur 
in Jfc^. It is given by 

V{r, A) = li^e-'^^ . (44) 

Let us recall that, according to the ETAS model, spontaneous and trig- 
gered events have statistically independent magnitudes whose PDF is to GR 
law 

^ _dQ(jnl ^ io-^("°)MnlO . (45) 
dm 

Thus, the spontaneous earthquakes have statistically independent magni- 
tudes {rrij}, j = 1, 2, . . . , r. According to the branching property of the ETAS 
model (Hawkes and Oakes, 1974; Helmstetter and Sornette, 2002), these 
spontaneous events trigger statistically independent numbers Rj(tk, t, r, mj, mo) 
of events within time interval [t,t + r]. This remark allows us to obtain the 
GPF Qki{z] T, mo) of the random number of events within [t, t + r] triggered 
by spontaneous events which occurred inside a given small interval 1^^: by 
replacing in ()43p the term by 

(46) 



(47) 



p(m')9(z; tk^ — t, r, m', mo)c?m' , 
using the expression ()44j) and performing the summation yields 



p{m')Q{z; tfci — t, T, m', mo)dm' 



mo 



fifci(^; T,mo) = exp i^uA 
Here, 

e{z; tk, , r, m', mo) = {z^^^''^ ,*>-,'»',mo)^ (4g) 

is the GPF of the random numbers of events triggered inside the time interval 
[t, t + r] by some spontaneous event of magnitude m' which occurred at a 
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time tk^ < t. Due to the statistical independence of the random numbers 
of events triggered by different spontaneous events which belong to different 
small intervals J^^, the GPF of the random numbers of events (including all 
their "aftershocks"), triggered within the interval + r] by spontaneous 
events which occurred before the window [t, t + r], is equal to 



Vti{z] T,mo) = '[[QkAz;tk^T,mo) = 
fcl 

exp uiA / p{m')<d{z;tk^ — t,T,m' ,mQ)dm' — 1 

\ I. Umn 



(49) 



In the limit A — > 0, the discrete sum in the exponential becomes an integral, 
which yields 



ill (z; r, mo) = exp ( uj 



dt' 



p{m')Q{z; t' — t, T, m', mQ)dm' 



mo 



(50) 

Analogously, one can find the GPF VL2{z\T,mQ) of the number of win- 
dowed events (including all their "aftershocks") triggered by spontaneous 
events which occurred within the interval + r]. It reads 



VL2{z] r, mo) = exp ( uj 
where 



t+T 



dt' 



p{m')zQ{z] t + T — t', m', mQ)dm' — 1 



mo 



(51) 



9(2;; r, m', mo) = 0(z; t = 0, r, m', mo) (52) 

is the GPF of the random number of events triggered inside the interval 
t G [0, r] by a spontaneous event which occurred at t = 0. 

For the third contribution from spontaneous sources in 1^^ G {t + T, 00), 
the causahty principle implies that = 1. 

Thus, the whole GPF of the numbers of events occurring within the time 
interval [t, t + r] is VL{z] t, mo) = ^^1(2;; r, mo)f22(^; t, "'^o); which can be writ- 
ten 



Q{z;t, mo) 



-LijL{z;T,mo) 



where 



L{z;T,mo) = / dm'p{m')x 

J mo 

/•oo rT 

/ [I — Q(z;t,T,m' ,mQ)] dt + / [I — zQ(z;t,Tn' ,TnQ)]dt 
Jo Jo 



(53) 



(54) 



In order to calculate GPF Q{z; ,r, mo) given by (jKHjl . (jSll), we need to 
determine the GPF G(z; t, r, m', mo) of the number of events triggered within 
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[t, t + r] by some spontaneous event magnitude m' which occurred at t' = 
0. Using a similar procedure as just done for fi(z;r, mo) and interpreting 
the Omori law as the PDF of random times of independently triggered 
aftershocks, we obtain 

Q{z;t,T,m' ,mo) = exp i — p{m') / dm"p{m")x 

\ J mo 

[b{t + r) - z<^{t + r) (g) e{z] r, m", mo) - $(t) ® Q{z; t, r, m", mo)] ^ . 

(55) 

The symbol denotes the convolution operator. We also define 

b{t) = 1 - a{t) , (56) 
where a(t) is given by (j2HI)- Similarly, Q{z; r, m', mo) satisfies to the equation 

Q{z; T, m', mo) = 

exp ( —p{m') I dm"p{m") [6(r) — z^{t) <d{z; r, m", mo)] ) . 

The above GPF's take into account events which have arbitrary magni- 
tudes m' ^ mo- In order to get the GPF which only counts observable 
events with magnitudes larger than some threshold level m, one has to re- 
place (i) mo by m in the r.h.s of expressions and and (ii) z by 

z{m' , m) = H{m — m') + zH{m' — m) , (58) 

where H{x) is unit step (Heaviside) function. We also replace the GPF 
0(2;; t, r, m', mo) by the GPF Q{z;t,T,m' ,m) of the number of observable 
aftershocks (whose magnitudes are larger than m), triggered by a mainshock 
of magnitude m'. We thus obtain 

n{z; r, m) = g-'^^^^'^'™) , (59) 



where 



L{z;T,m)= / [1 — D{z;t,T,m)]dt+ 
Jo 

Hi - D{z]t,m) + (1 - z)D+{z]t,m)]dt 
Jo 



and 



(60) 



D{z;t,T,m) = / Q{z;t,T,m',m)p{m')dm' , 

J rrn) 

D{z;T,m) = / Q{z;T,m',m)p{m')dm' , (61) 

J mo 

POO 

D^{z]T,m) = / <d{z;T,m',m)p{m')dm' . 

J m 
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It follows from and (jFTj) . after replacing the GPF 6(2;; , t, r, m', mo) 
by <d{z; ,t,T,m' ,m), the GPF 0(z;r, m',mo) by Q{z]T,m' ,171) and z by 
z{m",m) (|58|) . that the above auxiliary functions satisfy the equations 



D{z; t, r, m) = ^ 6(t + r) - (g) £'(2; t, r, m) 
+ r) ® r, m) + (1 - 2;)<l>(t + r) ® r, m) 



(62) 



and 



r, m) = \E' 6(r) - $(r) (g) D{z; r, m) + (1 - 2;)$(r) (g) D+{z; r, m) 



D+{z]T,ni) 
Here, 



6(r) - $(r) (g D{z; r, m) + (1 - 2;)$(r) (g D+{z; r, m) 



p{m')e^P^'^'^ydm' = Q(m)^[Q"^/^(m)i/] . 



(63) 



(64) 



For our following calculations, it is useful to expand "^{y) into a power series 
with respect to y: 

\E'(y) ^ 1 — ny + j3y^ — 77?/^ , (65) 



where 



n 



K7 



f3 = -\n 



7 , 



r(i-7) 



n^(7- 1)^ 
27(2 - 7) ■ 



(66) 



Recall that n is the so-called branching parameter (which we used to obtain 
expression (j211)), defined as the average of the total number of first-generation 
aftershocks triggered by some mainshock. 

Equations (jH^ and are rather complicated nonlinear integral equa- 
tions. For the rest of this section, we use them together with expressions 
()59|) . ()6U|) only to derive the average number of observable earthquakes. In 
the next sections, we exploit them fully to obtain the PDF of the inter-event 
times. 

As recalled for instance in the Appendix of (Saichev and Sornette, 2006a), 
the usefulness of the theory of GPF is to obtain simple expressions for the 
average number of observable earthquakes: 



(i?(t,r,m)) 



dVL{z; T, m) 



dz 



z=l 



(67) 
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Using dMl), (EDI), this leads to 

{R{t,T,m)) =ujQ{m)T + uj I M{t,T,m)dt + lu I M{T,m)dt , (68) 



where 



M{t,T, m) 



dD{z; t, r, m) 



dz 



M{t, m) 



dD{z; T, m) 



z=l 



dz 



z=l 



(69) 



The functions in ()69|) satisfy the following equations derived from ()62j) and 

M{t, T, m) - n^{t) (g) M{t, r, m) = 



n$(t + r) (g) M(r, m) + n [a(t) - a{t + r)] Q{m) 



and 



M(r, m) - n$(T) (g) M(r, m) = nb{T)Q{m) 
It follows from (ffn|l that the following equality is true 



(1 - ra) / M(t, r, m)(it = na(r) (g) M(r, m) + nQ{m)A{T) 
Jo 

Using this equality in (jUHj) . we obtain 

(i?(t,r,m)) = 



(70) 
(71) 

(72) 



Q(m) 



n 



1 — n 



g(m)A( 



1 — 77, 



a{T)®M{T,m)+ M{t,m)dt 



where A{t) is given by (jHOI)- 

Using equation (f7T|) . it is easy to show that 



(73) 



n 



1 — n 
This leads to 



a{T)®M{T,m)+ M{T,m)dt 



n 



1 — n 



Q{m)T 



n 



1 — n 



x> uj[m] 
1 — n 



uj{m) = ujQ{m) 



Q{m)A{T) . 

(74) 



(75) 



which has a physically intuitive interpretation: the average number of ob- 
servable events is equal to average number of observable spontaneous events 
uj{m)T within the interval of duration r, multiplied by the factor 1/(1 — n) 
taking into account all the triggered "aftershocks" of all generations. The 
result ()75|1 which can be obtained more directly (Helmstetter and Sornette, 
2002; 2003c) serves as a consistency check of our GPF formalism. 
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5 Statistics of recurrence times in the con- 
text of the ETAS model: accounting for 
Corral's empirical analyses 

We now use the formalism developed in the preceding section to derive 
the probability -P(t) given by (j^2|) that no event occurs with a time interval 
of duration r. Using (j^HI), (jEDl), (lEH), it is equal to 



where 



and 



P(r) = e-^^^"'™) , (76) 



L(r,m)=/ N{t,T,m)dt+ N^{t,m)dt (77) 
Jo Jo 



(78) 



N{t,T,m) = 1 - D{t,T,m) , 
N_{T,m) = 1 - D{T,m) + D^{T,m) . 

These functions have a clear probabilistic interpretation. For instance, N{t, r, m) 
is the probability that some spontaneous earthquake, which occurred at time 
t = 0, will trigger at least one observable event (which can be an aftershock 
of arbitrary generation) within the time window [t,t + r] {t > 0). 

It follows from that the probability N{t,T, m) satisfies the equation 

N{t,T,m) = 

(79) 

1 - ^ [$(t) (g) N{t, T, m) + <l>(t + r) (g) iV_(r, m)] . 
Analogously, it follows from (j63|l that N-{T,m) satisfies the equation 

A^_(r,m) = 

(80) 

1 — \1' [$ (r ) ® {j^m )] + Q(m)^ Q"^/^(m)$(r) ® N_ [T,mj\ . 

These are still rather complicated nonlinear integral equations. We can 
simplify them by linearization, which provides a reasonable approximation, 
as can be seen from the following argument. First, we notice that the in- 
equalities obviously hold 

N{t,T,m) < N{m) , N_{T,m) < N^{m) , (81) 

where 

N{m) = N{t = 0,T = oo,m) , N^{m) = N^{t = oo,m) . (82) 
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The majorants N{m) and A^_(m) satisfy the following nonlinear functional 
equations derived from from (f?^ . (jHUj) : 



N{m) = 1 - ^[A^_(m)] 
A^_(m) = 1 - $[A^_(m)] + g(m)^[g- 



/^(m)A^_(m)] . 



(83) 



We have solved numerically these equations (jH!^ and Fig. IH plots these solu- 
tions for N{m) and A^_(m). One can observe that, even for moderate mag- 
nitude thresholdss m — mo > 1, both N{m) <C 1 and N_{m) -C 1, so that, 
with the inequalities (jHT|) . one can be confident that a linear approximation 
of the r.h.s. of equations (fT9|l . (j8n|l should provide a good approximation. 

The linear approximation of equations (f7n|) and (fHUjl amounts to replace 
the function defined in ()fi4j) by the first two terms of the expansion 



This approximation ()84j) has a transparent probabilistic interpretation. It 
corresponds to assuming that any given earthquake can either trigger no 
or just one observable aftershock of first-generation. However, through the 
avalanche process of aftershocks which themselves trigger no or just one af- 
tershock, the total average number (R) = n/{l — n) of aftershocks over all 
generations which are triggered by a given mainshock can still be large for 
n close to 1. A possible justification of the approximation that no more 
than one observable first-generation aftershock is triggered by a given event 
is that this statement refers to observable aftershocks with magnitudes above 
the threshold level m > mg. For m large enough above mo, a given main- 
shock may trigger many aftershocks of first generation of magnitude larger 
than mo, but few of them will be observable with magnitudes larger than 
m, due to the fast decay of the GR PDF Q with increasing m. This effect 
has been shown to lead to a renormalization of the branching ratio n into 
an apparent value Ua < n (the equality holding only for the two fixed points 
n = and n = 1), characterizing the apparent branching structure of the 
triggering of observable events (Sornette and Werner, 2005b; Saichev and 
Sornette, 2006b). 

With the linear approximation ()84j) . equation ()8()j) takes the form 



^{y) ^ 1 



- ny . 



(84) 



N_{T,m) = Q{m) + 6^{t) (g) N_{T,m) , 



(85) 



where 



7-1 

6 = n[l — Q 'I (m)] . 
Similarly, expression ([7^1) allows us to obtain 

N{t, r, m) = n$(t) ® N{t, r, m) + n$(r + t) (g) A^-(r, m) . 



(86) 



(87) 
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Integrating the expression ()84|) term by term over time yields 

N{t, T,m)dt = n / N{t, r, m)dt + naij) ® (r, m) 



which simphfies into 

/ Nit, r, m)dt = air) ® (r, m) . (89) 

JO 1 — n 

Using ()89|1 . we rewrite the probabihty ()7(j|l in the form 

P(r)=exp a{T)®N(T,m)-uo N(t,m)dt], (90) 

V 1 — n JO / 

which now involves only the function N_{t,m). 

Thus, in order to find the explicit expression for P{t) from ()90p. we need 
to solve equation (jS^ . We can already state that its solution should obey 
the following limiting condition 

0(m) 

lim N_{T,m) = N^{m) = . (91) 



T— >00 



We now introduce the new function 

N (t m) 

9ir,m)= -1 ^ N4T,m) = N4m)[l~giT,m)]. (92) 

It is easy to show that g{T, m) satisfies the equation 

g{T, m) - 5^{t) (g) 5((r, m) = 5a(r) , (93) 
from which we derive the following relation 

5-1 r 

a(r) (g) g{T, m) = A{t) H — / g{t, m)dt . (94) 

JO 

Substituting in the second relation of (jHH) and using (jMj) yields 

/ N / uj(m)T uoim) , r , N , \ / N 

P(r) = exp ( -yry - g{t,m)dt\ , (95) 

where 

A = -i^--i- (96) 

1 — n 1 — 

and uj{m) = ujQ{m) is the average rate of observable spontaneous events. 
Recall that 5 is given by (jHF)|) . 
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A standard technique to obtain the solution of equation is to use the 
inverse Laplace transform. It gives 

g{T,m) = go (^5"^) , (97) 

where 

go{z, m) = -- r Im ( ^ ~ "^j^' ] e-^'dy (98) 

and 

^y,9) = y'eyTi-e,y) . (99) 

Due to the extremely slow decay of the function a(r) (j28|l for ^ 1, the exact 
solution (jHTjl of equation can be described accurately by the quasi-static 
approximation 

a(r)5 

Substituting this quasistatic approximation ()100|) in expression (|^ yields 

^, , / uj{m)T , ^ . r a(t)dt \ 

P(r) = axp [-AJ. _ ^„)A ^ j . (101) 

It is interesting to notice that, for 5 <^ 1, one may neglect 5 and replace the 
expression ()101|) by 

P(r) ^ exp ( -cj(m)r - n^^^A(r) ) , (102) 
\ 1- n 



which coincides with the simplified model probability (jHzj) derived in section 
3. The conditions under which this approximation holds can be visualized 
from Fig. El showing the dependence of of 6, taking as a function of m — mo 
for different values of 7 = b/a. This dependence results from the competition 
between the decay of the GR law controlled by its exponent b and the value 
of the productivity exponent a. 

Fig. El shows that the condition 5 <^ 1 needed to obtain ()102p requires 
that m — mo should not be too large, while Fig. jHhas shown that the linear 
approximation (jH^ used to obtain these results is all the more valid the larger 
is m—rriQ. It is thus worthwhile to come back to the more complete expression 
(jlOll) for P{t), which holds for arbitrary values of 6. It is convenient to 
rewrite (f7S|) so that the average number of windowed events is equal to 

X = A(m)r , A(m) = . (103) 

1 — n 
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With these notations, we have 



P(r) = Lp{x, m) = exp (^—r]x ^ ^ J div^ ^)dy^ , (104) 



where 



and 



^(->->^) = ^,. + (l_'5)(, + ,). • (105) 

ri = - -, i/= 1-nA, e = cAm =-^^. 106 

1 — 1 — n 

Correspondingly, the sought dimensionless PDF of the recurrence times de- 
fined by 

/(..m) = ;^^4¥^ (107) 



is given by 



/(x,m) = 

(108) 

(Ovil - 5)t~\x + e)^"^5(^(x, m, 6^) + [77 + //^((x, m, e)fj '^{x, m) . 

Similarly to the observation that the simplified version ()36|) for the PDF 
of recurrence times depends very slowly on the magnitude threshold m, the 
PDF of inter-event times ()108|) predicted by the ETAS model also depends 
weakly on m, if 6* <^ 1 and 7 is close to 1. Fig. IHl plots the dependence 
of the PDF ()1()8|1 as a function of the dimensionless inter-event time x for 
9 = 0.03, 7 = 1.2 and n = 0.9 for different magnitude thresholds m = 2, 4, 6. 
For such a span of magnitude thresholds, the average rate p03|) changes by 
many orders of magnitudes, while the corresponding PDF's are amazingly 
close each other. Fig. IHl indeed examplifies the very weak dependence of the 
PDF on m. 

Fig. 13 shows the empirical PDF's of the scaled inter-event times con- 
structed by Corral (2004a), as well as Corral's fitting curves using expression 
. We superimpose on these curves our theoretical prediction ()108|) , which 
actually provides a better fit to the data, especially for small x < 10~^, where 
our ETAS prediction correctly accounts for the impact of Omori's law, as al- 
ready mentioned in our discussion of the simplified model PDF (jHE))- We 
present the theoretical prediction p08|) for two different magnitude thresh- 
old levels m — rriQ = 2 and 6, showing the very weak dependence on the 
magnitude of completeness. These results suggest that the parameter 6 of 
the (bare) Omori law needs to be small (typically 9 < 0.1) to account for 
the data. Similar fits are obtained for the range of parameters 1 < 7 < 1.2 
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and 0.8 < n < 1, in agreement with bounds previously obtained from the 
prediction of the ETAS model on the distribution of seismic rates (Saichev 
and Sornette, 2006a). 

6 PDF of recurrence times for multiple re- 
gions: accounting for Bak et al. (2002) 's 
empirical analysis 

Until now, we have derived the prediction of the ETAS model for the 
statistics of inter-event times in a single homogeneous region and have shown 
that it is compatible with the empirical observation of an approximate uni- 
fied scaling law of the form We have also shown that our theoretical 
expression fits remarkably well the empirical PDF's over the whole range of 
recurrence times, accounting for different regimes by using only the physics 
of triggering quantified by Omori's law. 

There is another unified saling law for the statistics of recurrence time 
obtained after averaging over multiple regions, which has been proposed by 
Kossobokov and Mazhkenov (1988) and Bak et al. (2002). We will refer for 
short to this law as Bak et al.'s unified law. The aim of this section is to 
derive analytically Bak et al.'s unified law obtained for multiple regions based 
on the prediction of the ETAS model that we have presented in the previous 
section and on the simple hypothesis that all regions obey the PDF form 
()l()8j) but with different average seismic rates. This simply means that Bak 
et al.'s unified law in our view just results from an appropriate averaging 
of expression ()108|) over the statistics of the average rates of the multiple 
regions under observation. 

We first present in subsection 16.11 a general scheme to perform this aver- 
aging in the next subsection with various test statistics for the average rates 
in different regions. Then, in the second subsection 16. 2[ we introduce a frac- 
tal geometrical model which allows us to propose a specific statistics for the 
average rates in different regions, which lead to excellent fits to the empirical 
data. 

6.1 General statistics for the average regional seismic 
rates 

In order to choose reasonable statistics for the average seismic rates cjj 
in multiple regions with the minimum of additional parametes, we assume 
that these uji^s obey a principle of statistical self-similarity: if the branching 
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ratios rij's and the linear size Lj's of different regions are identical, then the 
seismic rates can be written 



u}i = uj Ui , (f09) 

where uj is the rate of spontaneous events averaged over all regions, while the 
Mj's are mutually independent random variables distributed according to the 
same PDF £{u). The PDF £{u) is assumed to be independent of the common 
values of the branching ratio n and linear size L of the multiple regions. 
Additionally, we assume that, as for the scahng law © of a single region, the 
PDF £{u) does not depend on the magnitude threshold m (which is common 
to all regions). One may interpret this assumption as a consequence of the 
scaling properties of the GR law. 

Since the average rates of spontaneous events cannot be directly observed, 
it may be preferable to replace the statement by the equivalent repre- 

sentation 

Ai = A Mi , (110) 

which involves the total observable seismic rates above a magnitude thresh- 
old. Due to pi()|l . the PDF £{u) should satisfy to double normalization 
condition 

/ £{u)du = / u£{u)du = 1 . (Ill) 
Jo Jo 

In this subsection, we concentrate our attention on the universal properties 
of the PDF of recurrence times constructed over multiple regions, which 
are found almost the same for qualitatively different distributions £{u). In 
the next subsection l(j.2| we will obtain a specific expression for the PDF 
£{u) from a simple model based on a fractal distribution of spontaneous 
earthquake sources. 

Applying our conjecture ()11U|) to the general relation ()1U|) . the PDF of 
inter-event times over multiple regions is obtained as 

It is convenient to represent the probability -P(t) given by (fTBj) that no events 
occur in a given single region in the form 

P(r) = e-"^(^) , (113) 

where the auxiliary function 

K{x) =iuL(^j,m}j , x = Xt, (114) 
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does not depend on u. It follows from (jlOU) that 

rx 

K{x) = rix + iy g{y,m,6)dy , (US) 
Jo 

where g{y,m,9) is given by Substituting into ()112j) yields the 

following expression analogous to for the PDF of the recurrence times 
over multiple regions: 

H(r) = XhiXr) (116) 

where 

We have used the following Laplace transforms 

£{x) = / £{u)e-''''du , £i{x) = ^ = / u£{u)e-''''du . (118) 

JO dx Jo 

It is possible to derive the following asymptotic behavior of the PDF h{x) 
()117p . Below, we will illustrate the main features of these asymptotics using 
the limiting case 6 = and r 3> c such that K{x) reduces to 

K{x) ^ {l-n)x + j:e^x^-\ (119) 

1 — 

corresponding to the simplified model of recurrence time statistics. 

Let us first consider the behavior of h{x) for x -C 1. First, the asymptotic 
behavior of the dimensionless PDF h{x) for x <^ 1 is universal in the sense 
that it does not depend on the shape of the PDF S{u) of the average rates 
of the multiple regions. Indeed, due to the double normalization condition 
(jllH) . we have 

£{K) ~ £i{K) ~ 1 (fsT < 1) . (120) 
Thus, in view of ()117j) and (jll9j) . we obtain 

/i(x) ^ (1 -n + ne V^)^ + ^ne V^-^ (x < 1) . (121) 

For very small x's, the last term prevails and we have the universal asymptotic 

h{x) ^ ene^x~^'\ x < < 1 . (122) 

The range of validity of this asymptotic can be roughly estimated by putting 
formally e^x~^ ^ 1 in the r.h.s. of relation p2H) . which gives 

h{x) ^ 1 + 9nx-^ (123) 
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so that 

x^^On . (124) 

Let us now consider the behavior of h{x) for x ^ 1. We observe an 
almost universal asymptotic of the dimensionless PDF h{x) given by pi7|l 
for X ^ 1. Indeed, let us assume that the PDF of the average seismic rates 
of the multiple regions has, for small u, the following power asymptotic 

£{u) ^ /iu" , u < 1 . (125) 

This leads to the explicit form of the Laplace transforms (jll8|l : 

^(a;) ~ /ir(l + a)x~^"" , ^i(x) ~ /ir(2 + a)a;"2~" , x > 1 . (126) 

This finally yields 

h{x) ^ ^r{l + a)x 

(127) 

;i + a)(l - n + ne^x-^K-'^-'^ix) + One^x-^-^R-^-'^ix) . 

For small 9, the leading behavior of this asymptotic can be obtained by using 
K{x) ~ X, which gives 

h{x) ~ . (128) 

A useful overview of the behavior of the dimensionless PDF h{x) valid for 
all interdiate values x, including the two above asymptotics p22|l and p28j) . 
can be obtained by considering a smoothed function interpolating between 
them. Let us for instance consider the following modeling function 

, , , 1 + Onx'^-^ 

hm{x) = -Y^^^ . (129) 

We put here for simplicity fiT{2 + a) ~ 1. Fig. |H1 shows the function xhm{x) 
for 6 = 0.03, n = 0.9 and for different values a, which mimics the empirical 
construction reported by Bak et al. (2002). Note in particular the existence 
of a plateau for x <^ 1, corresponding to the Omori law asymptotic p22|) . 
One can also observe the power law asymptotic corresponding to the 

asymptotic ()128|1 for x ^ 1. In addition, there is an intermediate kink for 
X ~ 1, due to the first term of the asymptotic formula ()123|1 . All three 
features mimic the behavior found by Bak et al. (2002). 

These properties of h{x) are weakly dependent on the value of the ex- 
ponent a, which controls the asymptotic behavior p25|) of the PDF of the 
average seismic rates of the multiple regions for small rates. To illustrate this 
point, consider two qualitatively different sample distributions S{u) which 
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obey the double normalization Our first example is the Gamma dis- 

tribution formulated so that it obeys 

n _|_ rv)^^" 

SS{u) = V.T N ^"e-(^+")" . (130) 
r(l + a) 

The main properties of this distribution are the presence of a power asymp- 
totic ()125|) for small u's and an exponential decay for u ^ 1. In this case, 
the Laplace transforms ()118|) are 

^ ^ Vl + a + x/ ^ \l + a + xj ^ ' 

Our second example for a distribution obeying ()111|) is 

£» = P±ifl + !iV, (132) 



which corresponds to the particular case a = 0. The main difference between 
this distribution and the previous Gamma distribution is the slowly decaying 
power law tail for m ^ 1: 

E^[u) ~ (m > 1) . (133) 

In this case, the Laplace transforms piSj) read 

^P(x) = (1 +p)(px)^+^'e*'^r(-l -p,px) (134) 

and 

El{x) = -(l+p-(l+p + px)£:„(x)) . (135) 

X 

Fig.Elplots the PDF h{x) obtained with the two examples ()13()|1 and ()132|1 . It 
would be difficult to distinguish between these two distributions in a realistic 
empirical setting with limited data and large statistical fiuctuations. 

In the next subsection 16.21 we develop a simple model for E{u) based 
on the idea that seismic rates of spontaneous events are distributed on a 
fractal geometrical set. In a nutshell, the key idea is the following: a spatial 
distribution of spontaneous events with fractal dimension d <2 implies that, 
the larger is the number of events, the larger are the spatial distances between 
them. Thus, the PDF of the number R of events within a given region 
should decay quite fast at the number of events increases. This tends to 
favor the decay for large seismic rates modelled by the Gamma distribution 
p3()j) over the power law distribution p32j) . To finish this general discussion. 
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we present in Fig. the empirical data analyses by Corral (2004a) together 
with prediction for h{x) obtained by using the Gamma distribution for S{u) 
with 6 = 0.03, n = 0.9, e = 0.76 and a = 0.2. The fits to the empirical data 
are suggestive as they captures all its qualitative features. There are however 
undeniable quantitative differences which result in part from the fact that 
the Gamma distribution is probably not the true distribution for S{u) even 
if it captures correctly the mean properties of the data. 

6.2 Probabilistic consequences of the fractal geometry 
of the spatial distribution of spontaneous earth- 
quakes 

In the framework of the ETAS model combined with the hypothesis of 
statistical self-similarity, the preceding subsection 16.11 has stressed that an 
analytical understanding of Bak et al.'s unified law requires the knowledge 
of the PDF S{u) of the average rates over the multiple regions. To constrain 
its structure, we make the assumption that the well-known fractal spatial 
organization of all observed earthquakes reflects a similar fractal spatial or- 
ganization of the subset of spontaneous events which are the origin of the 
earthquake triggering activity. We thus assume that spontaneous events are 
distributed spatially on a fractal geometry, which we are going to exploit to 
derive a natural form for the PDF S{u). 

For this goal, the key idea is based on an important relation between two 
reciprocal random variables. The first random variable is the area S{k) occu- 
pied by k spontaneous events which occurred during a fixed time interval of 
duration r. The second random variable k{S) is the inverse function of S{k), 
which is nothing but the number of spontaneous events within the given area 
5*. Although both functions are ill-defined in practice due to the difficulties 
in defining unambiguously the areas occupied by chaotically occurring dis- 
seminated events, we can nevertheless obtain approximate scaling laws for 
the distributions of these two variables, which derive from the fact that the 
average k{S), which is proportional to the rate of spontaneous events uj{s), 
should obey to power law 

{k{S)) UJ{S) r-. L'^ . (136) 

Here, L is a characteristic scale of the area S and d is the fractal dimension 
of the spatial set of spontaneous events. 

In order to obtain these scaling laws, we recognize that the i-th sponta- 
neous event "occupies" a spatial area Si such that, if some area 5* contains k 
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spontaneous events, then it is equal to 

S{k)^s, + S2 + --- + Sk . (137) 

One may interpret Si as given by Sj ^ Lf where Li is the distance between the 
2-th event and its nearest neighbor. Next, let us assume that the summands 
in ()137p are statistically independent with a common PDF v{s) which has 
the following aymptotic power law 

vis) ^ ' ''^'^ 0<9<l. (138) 

The expression p38|l corresponds to a distribution of spatial jumps, as in a 
so-called Levy flights (Metzler and Klafter, 2000), which underlies the fractal 
spatial distribution of spontaneous events. From the generalized central limit 
theorem (Gnedenko and Kolmogorov, 1954; Sornette, 2004), the PDF of the 
random sum ()137|) for A; ^ 1 is given asymptotically by 

where fg{x) is an infinitely divisible distribution, whose Laplace transform is 
equal to 

faiy) = / fg{x)e-y^dx = e-^' , 0<g<l. (140) 

Let us now view expression ()137|1 from a different vantage. One can also 
interpret ()137|1 as an implicit equation 

S{k) = S (141) 

for the unknown variable equal to the number k{S) of spontaneous events 
occurring within the given area S. Piryatinska et al. (2005) have shown 
that the PDF of the random variable S allows one to derive the PDF of the 
random variable k{S) as 



WftS) = |;uJ|| , (142) 



where 



gxW/g-'y y^i/g 



We then deduce that the average of the random variable k{S) distributed 
according to the PDF p42|l is given by 
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Since ~ in p44|) and comparing with ()136|) yields the self-consistent 
determination of the exponent g 



g = d/2 



(145) 



which is indeed between and 1 for < d < 2. Assuming that average rates 
of single regions are proportional to k{S), that is A ~ c<j ~ k{S), we obtain 



Ui 



iHS)) 



(146) 



where Si are different regions with the same linear size L. The sought PDF 
of the average seismic rates is thus equal to 



S{u) = -Wd/2 (- 

u \u 



u = T{l + d/2) . 



(147) 



One can show that 



Wg{u) 



1 

Tcg Jo 



exp 



x^^^ cos 



7r 



cos 



ux — x^^^ sin 



dx 



0.5<g <1 . 

(148) 

Fig. ^2 plots the PDF S{u) given by ()148j) for different values of the fractal 
dimension d of the set of earthquake epicenters. The tail of S{u) is more 
extended for smaller fractal dimensions d. 

The explicit form p47|) of S{u) allows us to obtain the Laplace transforms 
S{x) defined by pi8|) as 



S(x) 



E, 



d/2[ 



-UX] 



(149) 



where Eg{z) is the Mittag-Leffler function which has the following integral 
representation 



TC 



sm{ng) 



^dy 



v'^ + y'^s -\- 2vy3 cos{'Kg) 



(150) 



We then calculate explicitely the corresponding PDF h{x) given by ()117|) 
of the inter-event times associated with ()147|1 . The results are shown in 
Fig. El for three different values d = 1.2, 1.4, 1.6 of the fractal dimension of 
earthquake epicenters. 

Fig. El is the culmination of the present work. It first makes use of 
our theoretical analysis based on the ETAS model leading to the prediction 
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(jlOSp for a single homogeneous seismc region, which is combined for mul- 
tiple regions with different seismic rates distributed according to the PDF 
()147|1 predicted from our fractal model of earthquake epicenters, to finally 
obtain the global PDF ()117|1 of inter-event times. Our final prediction ()117|) 
is compared to the empirical distribution represented in Fig. 2 of (Corral, 
2004a), for the parameters d = 1.8, n = 0.9, 7 = 1, 1 and to 3 values of the 
magnitude threshold m — tjiq = 2; 4; 6. The agreement is remarkable. 

7 Concluding remarks 

We have proposed a general theoretical approach to describe the statis- 
tics of recurrence times between successive earthquakes in a given region, or 
averaged over multiple regions. This work was motivated by the reports by 
several authors that recurrence times between successive earthquakes should 
be considered for broad areas, rather than for individual faults, and could pro- 
vide important insights in the physical mechanisms of earthquakes. Several 
authors mainly from the Physics literature have proposed that the scaling law 

found to describe well empirical data reveals a complex spatio-temporal 
organization of seismicity, which can be viewed as an intermittent flow of 
energy released within a self-organized (perhaps critical) system, for which 
concepts and tools from the theory of critical phenomena can be applied. 

We have shown that this view is probably too romantic because much 
simpler explanations can be proposed to fully account for the empirical ob- 
servations. Indeed, we have shown that the so-called universal scaling laws of 
inter-event times do not reveal more information than what is already cap- 
tured by the well-known laws of seismicity (Gutenberg-Richter and Omori, 
essentially), together with the assumption that all earthquakes are similar 
(no distinction between foreshocks, mainshocks and aftershocks). This con- 
clusion is reached by a combination of analyses, which start from a gener- 
alization of Molchan (2005) 's argument, to go to simple models of triggered 
seismicity taking into account Omori's law, and end with a detailled study of 
what the ETAS model has to tell us on the statistics of inter-event times. By 
using the formalism of generating probability functions, we have been able 
to derive analytically specific predictions for the PDF of recurrence times 
between earthquakes in a single homogeneous region as well as for multi- 
ple regions. Our theory has been found to account quantitatively precisely 
for the empirical power laws found by Bak et al. (2002) and Corral (2003; 
2004a). We showed in particular that the empirical statistics of inter-event 
times result from subtle cross-overs rather than being genuine asymptotic 
scaling laws. We also showed that universality does not strictly hold. 
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Therefore, to the question raised in the introduction, we are led to con- 
clude that the statistics on inter-event times described in (Bak et al., 2002; 
Corral, 2003, 2004a,b, 2005a; Livina et al, 2005) is not really new in the 
sense that they do not reveal information which is not already contained in 
the known laws of seismicity. Our conclusion is that they can be derived from 
the known statistical properties of seismicity, so that they are only different 
ways of presenting the same information. In particular, the fact that the 
simple models of Lindman et al. (2005) are not able to fully reproduce the 
structure of the empirical statistics of recurrence times, as noted by Corral 
and Christensen (2006), does not necessarily imply that the so-called "uni- 
fied scaling laws" reveal any novel information. We think we have clearly 
shown that they can be derived from simple and well-known physical ingre- 
dients, when carefully taking into account the physics of triggering between 
earthquakes. In this sense, the present work is the continuation of an effort 
to classify the empirical observations which can be, from those which cannot 
be, explained by the simple ETAS benchmark (Helmstetter and Sornette, 
2003a; 2003b; Saichev and Sornette, 2005; 2006a). With this effort, we hope 
to eventually help identify real robust statistics which can not be explained 
by the ETAS benchmark or its siblings, leading us towards the acquisation 
of new interesting and important knowledge on the physics of earthquakes 
with potential applications for their forecasts. 
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Figure 1: Plots of the PDF f{x) §^ ioi n = 0.9, 9 = 0.03, for different 
tfireshold magnitudes m — mi = 0; 2; 4; 6 (top to bottom on the left), where 
the reference parameter is ei = Aic = 10~^ for mi. 
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Figure 2: Plot of the PDF f{x) given by expression (jHUj) for the parameters 
n = 0.9, 9 = 0.03, ei = 0.76, m — mi = 2, and Corral's fitting curve (jHHj) 
for 7 = 0.38, 5 = 1, d = 1.7 and C = 0.75. We also show an intermediate 
asymptotics ~ discussed by Corral (see below) as well as the short- 

time asymptotics ~ corresponding to the Omori law. 
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Figure 3: Partition of the time axis in small intervals which are classified in 
three categories Ik^ , Ifc^ and Ik-^ . 
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Figure 4: Plots of the solutions of the nonlinear functional equations ()83|) for 
N{m) and N_{m), for the parameters 7 = 1.2 and two values of criticality 
parameter n = 0.8; 0.9. 
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Figure 5: Dependence on m of the parameter 6 defined in (jH^ . which resuhs 
from the competition between the GR and productivity laws. Bottom to top: 
7 = 1.1; 1.2; 1.3. 
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Figure 7: Plots of the empirical PDF's of scaled inter-event times taken 
from (Corral, 2004a), Corral's fitting curves using expression ()H8j) (single solid 
curves) and our theoretical prediction ()108|1 for = 0.03, 7 = 1.2 and n = 0.9 
(double solid curves). Our ETAS prediction is shown as double solid curves 
because we show the theoretical prediction ()108p for two different magnitude 
threshold levels m — mo = 2 and 6, showing the very weak dependence 
on the magnitude of completeness. Top to bottom: the NEIC worldwide 
catalog for regions with L > 180 degrees, 1973-2002; NEIC with L < 90 
degrees, (same period of time); Southern California, 1984-2001, 1988-1991, 
and 1995-1998; Northern California, 1998-2002; Japan, 1995-1998, and New 
Zealand, 1996-2001; (bottom), Spain, 1993-1997, New Madrid, 1975-2002, 
and Great Britain, 1991-2001. The curves are translated for clarity by the 
factors 1; lO'^; lO^^; iQ-^; iQ-s and 10"^° from top to bottom. 
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Figure 8: Plots of the function xhm{x) for 6 = 0.03, n = 0.9 and for different 
values of a = 0; 0.1, 0.2, demonstrating the approximate unified scaling law 
of the PDF of inter-event times obtained over multiple regions, as reported 
by Bak et al. (2002). 
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Figure 9: Plots of the PDF h{x) given bv dTTTI) for 6 = 0.03, n = 0.9, obtained 
by using the two examples ()130p and p32|) for the the PDF of the average 
seismic rates of the multiple regions. Solid line: the Gamma distribution 
()130j) for a = 0; dashed line: the power tail distribution ()132|1 for p = 0.5. 
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Figure 10: Empirical PDF's of the recurrence times between earthquakes 
over multiple regions following Bak et al. (2002) 's procedure obtained from 
Fig. 2 of (Corral, 2004a), on which has been superimposed our prediction for 
h{x) obtained with the Gamma distribution ()130|1 for S{u) with 6 = 0.03, 
n = 0.9, e = 0.76 and a = 0.2. The curves have been translated from top to 
bottom by the factors 1; 10"^; 10"^ 10"^ 
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Figure 12: Plots of the PDF h{x) of inter-event times given by ()117|1 cor- 
responding to the PDF p47|l of the average regional seismic rates. Top to 
bottom on the right side of the picture d = 1.2; 1.4; 1.6. The straight lines 
show asymptotic power laws, as proposed by Bak et al. (2002). 



50 



jQ-r> JQ-S ^Q-l fQ-a JQ-2 IQ-l 1 IQ J(j3 

Figure 13: Plots of the empirical PDF of inter-event times in multiple regions 
collected in Fig. 2 of (Corral, 2004a), which is compared with our theoret- 
ical PDF h{x) given by ()117|1 (thick lines), represented for the parameters 
d = 1.8, n = 0.9, 7 = 1,1 and three values of the threshold magnitude 
m — mo = 2; 4; 6. Corral calculated the recurrence-time probability densi- 
ties following Bak et al. (2002) 's procedure, after rescaling by the regional 
seismic rates. 84 distributions are shown, with region sizes L ranging from 
0.039° to 45°, and threshold magnitudes m varying between 1.5 and 6. The 
curves are shifted by factors 1; 10~^; 10~^ and 10~^, and correspond to: 1 
(top). Southern California, 1984-2001; 2, Northern California, 1985-2003; 3, 
Southern California, 1988-1991 (stationary rate), NEIC, 1973-2002, Japan, 
1995-1998, and Spain, 1993-1997; 4 (bottom). New Zealand, 1996-2001, and 
New Madrid, 1975-2002. The thin straight lines are Corral's fits by power 
laws of the short or intermediate times: 1 (top), 0.12/x'''^^; 2, 0.15/x°'^, 3 
and 4 (bottom), 0.05/x°'^^ and 0.5/x°'^. In all cases. Corral has fitted the 
long-time tail by 0.25/x^'^. The times in the horizontal axis span from 2 min 
to about 20 years. Recurrence times smaller than 4, 10, and 2 min are not 
shown, for Japan, NEIC, and the rest of catalogs, respectively. 
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